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Loopy belief propagation (LBP), which is equivalent to the Bethe approximation in statistical 
mechanics, is a message-passing-type inference method that is widely used to analyze systems based 
on Markov random fields (MRFs). In this paper, we propose a message-passing-type method to 
analytically evaluate the quenched average of LBP in random fields by using the replica cluster 
variation method. The proposed analytical method is applicable to general pair-wise MRFs with 
random fields whose distributions differ from each other and can give the quenched averages of 
the Bethe free energies over random fields, which are consistent with numerical results. The order 
of its computational cost is equivalent to that of standard LBP. In the latter part of this paper, 
we describe the application of the proposed method to Bayesian image restoration, in which we 
observed that our theoretical results are in good agreement with the numerical results for natural 
images. 
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I. INTRODUCTION 

Loopy belief propagation (LBP) [I|, which is a 
message-passing-type inference method, is widely preva¬ 
lent in various fields, including computer science, as a 
powerful tool for statistical procedures in systems based 
on Markov random fields (MRFs) ii- LBP is equiv¬ 
alent to the Bethe approximation in statistical mechan¬ 
ics iS and is also known as the cavity method. An 
analysis of the statistical behaviors of LBP is important 
to develop an understanding of LBP. In this paper, we 
focus on LBP in pair-wise MRFs with random fields and 
we present a statistical analysis of it, namely, an analysis 
of the quenched average of LBP over random fields. The 
topic of pair-wise MRFs in random fields is an important 
research field in statistical mechanics As described 

below, a statistical analysis of LBP in random fields is 
also important for the field of Bayesian signal processing 
in computer science. 

Bayesian image restoration Q, in which images de¬ 
graded by noise are restored using the Bayesian frame¬ 
work, is an important generic technique for various types 
of signal processing. Suppose that there is an original 
image and that the original image is degraded through a 
specific noise process. We observe only the degraded im¬ 
age as the input, and we want to produce the restored im¬ 
age as the output. From the statistical mechanics point of 
view, the standard framework of Bayesian image restora¬ 
tion corresponds to the framework of a two-dimensional 
ferromagnetic spin model in random fields BIS- In this 
correspondence, the input image, namely, the degraded 
image, is regarded as the random fields in the Bayesian 
image restoration system. 


Since the model used in the Bayesian image restora¬ 
tion system is designed by using an intractable pair-wise 
MRF, LBP is often applied to implement it. Hence, in 
the evaluation of the statistical performance of the im¬ 
plemented image restoration system, we encounter the 
evaluation of the quenched average of LBP over the ran¬ 
dom fields, namely, over the input images. For this pur¬ 
pose, for Ising systems, the authors proposed an ana¬ 
lytical evaluation method for it . In the previous 

method, the evaluation of the quenched average of LBP is 
reduced to solving simultaneous integral equations with 
respect to the distributions of the messages. However, 
the method is not very practical, because the computa¬ 
tional cost of solving the integral equations is consider¬ 
able and its approximation accuracy is poor. Further¬ 
more, the method cannot evaluate the quenched average 
of the free energy and is formulated only in Ising systems. 

In this paper, we propose a new analytical method for 
evaluating the quenched average of LBP over random 
fields based on the idea of the replica cluster variation 
method (RCVM) The presented method allows 

the quenched average of the Bethe free energies over ran¬ 
dom fields in general pair-wise MRFs to be evaluated, 
unlike the previous method. 

The remaining part of this paper is organized as fol¬ 
lows. A brief explanation of LBP is given in section HIl 
Section imi constitutes the main part of the paper. The 
proposed method is shown in section IHI Ci and some 
numerical results for checking its validity are shown in 
section UlI PI In section UlI E[ we show a case that is ex¬ 
actly solvable by the present method. In section ITVl we 
explain the framework of the framework of Bayesian im¬ 
age restoration and compute the statistical performance 
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of the Bayesian image restoration system using the pro¬ 
posed method. Finally, section |V] closes the paper with 
concluding remarks. 


and inverse temperature, which takes a positive value, re¬ 
spectively. In this paper, only random fields h are treated 
as the quenched parameters. 


B. Loopy Belief Propagation 


II. LOOPY BELIEF PROPAGATION IN 
RANDOM FIELDS 

A. Model Definition 

Let us consider an undirected graph G{V,E) consist¬ 
ing of n vertices and some edges, where V = {1, 2,..., n} 
is the set of vertices, and E = {{i,j}} is the set of edges 
between a pair of vertices, where {i,j} denotes the undi¬ 
rected edge between vertices i and j. On the graph, with 
the discrete random variables S' G {S'i | i e E}, let us 
define the pair-wise MRF expressed by 


It is known that LBP is derived from the minimum 
condition of the variational Bethe free energy Q. In 
this section, we give a brief explanation of the deriva¬ 
tion of LBP according to the cluster variation method 
(cvM) [il,[i^. The free energy of the MRF in equation 
(ED is defined by 

F(/i,/3) :=J2H{S-h)P{S\h,f3) 

s 

+ l^P(S|/i,/3)lnP(S|/i,/3). (2) 

^ s 

In the Bethe approximation in the CVM, we approximate 
the MRF by 


P{S\h,f3):=j^exp{-(iH{S;h)), ( 1 ) 

where 

i€V 

is the Hamiltonian of the MRF. Here, hi) is a spe¬ 

cific function of the variable Si and the random field hi 
on vertex i, and (j)ij{Si, Sj) is a specific function on edge 
{i,j}- The notations, Z and /3, are the partition function 


P{S\h,f3) 


(Oiey bi{Si)) ( n{jj}e£; bi,j{Si, Sj)) 


( 3 ) 


where bi{Si) and bij{Si, Sj) are the one-vertex and two- 
vertex marginal distributions (or the beliefs) of the MRF. 
This approximation corresponds to the cluster decompo¬ 
sition shown in figure [T] The right-hand side of equation 
(I3|) is the product of the marginal distributions of the 
clusters divided by the product of the marginal distribu¬ 
tions of the double-counted clusters. By applying this 
approximation to P{S \ h) in the logarithmic function of 
the last term in equation ([2|), we obtain the variational 
Bethe free energy expressed by 


x- !tu h 11_ iTf a. a i h ( HieV *^f)) 

.^bethe[{^i; bi,j}] ■— / . h)P[S \ h, f3) + / ^ P[S \ h, f3)\n h f Q \h f Q \ 

= —'^^'^^<f>i{Si,hi)bi{Si)— ^j)bi,j{Si, Sj) + — 'y^,'Hi[bi] 

iev Si {i,j}eESi,Sj ^ i&V 

{'^2[bi,j] - Uilbi] - 'Hi[bj]), 




( 4 ) 


IhllGB 


where 


and 


ni[bi\ ■.^Y.^i{Si)\nbi{Si) 

Si 


'^2[bi,j] ■■= bij{Si,Sj)\nbij{Si,Sj) 

Si,Si 


are the one-vertex and two-vertex negative entropies. 

LBP is obtained by the variational minimization of 
the variational Bethe free energy with respect to the 
beliefs. By minimizing the variational Bethe free en¬ 
ergy under the normalizing constraints, Es bi{Si) = 
Ssi Sj the marginalizing con- 
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(a) (b) 


FIG. 1. (a) Square graph with four vertices, (b) Cluster 

decomposition of the Bethe approximation of the CVM for 
(a). The vertices and two connected vertices are selected as 
the clusters. 

straints, bi{Si) = 'S'j) and bj{Sj) = 

Xs bij{Si, Sj), we obtain the message-passing equation 
of LBP: 

Mi^j{Sj) oc '^expf3{(j)i{Si,hi) + Sj)) 

Si 

X Yl (5) 

k£di\{j} 

where Mi^j{Si) is the message (or the effective field) 
from vertex i to vertex j. Using the messages satisfying 
the message-passing equation, we can compute the beliefs 
that minimize the variational Bethe free energy as 

bi{Si) (xexjp(l3(j)i{Si,hi)) Y\_ ( 6 ) 

jedi 

bij{Si,Sj) oc expP[4’i{Si,hi) + (j)j{Sj,hj) 

+ Sj)) 

kedi\{j} ledj\{i} 

(7) 

where di is the set of vertices connected to vertex i: di = 
{j I J S U, {i,j} € E}. The Bethe free energy of the 
MRF is the minimum of the variational Bethe free energy, 

Sbethe(^;/^) •— min .Fbethe[{bj, 62 ,j}]; 

{bi,bij} 

and is obtained by substituting the beliefs obtained by 
equations m and © into the variational Bethe free en¬ 
ergy in equation (U). In LBP, the beliefs obtained by 
equations ( 0 ) and © are regarded as the Bethe approx¬ 
imations of the true marginal distributions of the MRF. 
When an undirected graph G{V,E), on which the MRF 
is defined, has no loops, the Bethe free energy and the 
beliefs are equivalent to the true free energy and the true 
marginal distributions of the MRF, respectively. 

The main proposal presented in this paper is a method 
for evaluating the quenched average of the Bethe free 
energy over the random fields; 

[Ehetheih, P)]h := j dh(^Y[Piihi)^Ebetheih,l3), (8) 

•7 ^ iev 


where the notation represents the average value 

over the random fields, and Pi{hi) is the distribution of 
the field on vertex i, where these distributions can vary 
by vertex in general. 

III. PROPOSED METHOD 

According to equation ([5]), in principle, we have to 
perform the averaging operation after constructing the 
Bethe free energy to obtain [Fbethe]^ (see figure UKa)). 
However, it is not straightforward to directly integrate 
the Bethe free energy. Thus, we adopt another strategy. 

In this paper, we propose an approximate method 
based on the idea of the RCVM [l^ [l^. Figure [5Kb) 
shows the procedure of our method. In the method, we 



(a) (b) 


FIG. 2. Procedures of obtaining the quenched average of the 
Bethe free energy, (a) Principled procedure, (b) Procedure 
of the proposed method based on the RCVM. 

first take the average of the free energy using the replica 
method and CVM, namely, the RCVM fsection IlII Al) . 
and then, we apply the Bethe approximation to the re¬ 
sulting form of the RCVM fsection rill Bl) . If the exchange 
of the order of the two operations, the Bethe approxima¬ 
tion and the quenched averaging operation, is allowed, 
we can expect to obtain a good approximation of the 
quenched average of the Bethe free energy. 

A. Replica Cluster Variation Method 

First, we obtain the quenched average of the true free 
energy of the MRF in equation ©, that is, 

[F{h,l3)]h = -^ j dh(^l[p,ih,)) \nZ{h,f3). 

In the context of the replica method 1,113 , we have 

[F{h,P)]h = (9) 

P X 

[ dh(^l[Mhi))z{h,pr. 

•7 iev 


where 
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By assuming that x is a natural number, we obtain 

X 

Z,=^exp^(^e,(5,)+ ^ ^^-..,(^“,5“)), 

Sx iGV {iJjGE a—1 


is the Gibbs distribution of the x-replicated system, and 
5“ = {S'" I i € V}. The energy function HintiS) is 
dehned by 


where = {S'" | i e F, a = 1, 2,..., x}, = {S'" | a = 

1 , 2 ,..., x|, and the function eASi) is defined by 

et{Si) := gin / dh p^h) exp (^13 Y ^ ■ bjfe-® 


We regard as the partition of the x-replicated system 
and define the x-replicated free energy as 


Fx ■■= --InZx 

X 

= - E E + E E ifint (5“)Q“(S'“) 


lev Si 
1 


a=l S“ 


-a Y ^^Px{Sx), 


( 10 ) 


where 


1 ^ 

PxiSx) := — exp/3(Eez(5.)-E^-t(‘^' 


iev 


and is just the interaction term of the original system. 
The distributions, Qi{Si) and Q‘^{S°‘), are the marginal 
distributions of the distribution Px{Sx)- The factor 
graph representation of the x-replicated system is shown 
in figure EDa). 


In accordance with the cluster decomposition based 
on the CVM shown in hgureEKb), by using the marginal 
distributions, {Qi{Si),Q°'{S°‘)}, together with the one- 
variable marginal distributions of PxiSx), {QfC'S'f)}, we 
approximate the Gibbs distribution of the x-replicated 
system as 


Px{Sx 


(n.gy ni=i Qfjsn) (n.gy ( ni= 


As in equation ([3]), PxiSx) is approximated by the prod¬ 
uct of the marginal distributions of the clusters divided 
by the product of the marginal distributions of the 
double-counted clusters. By applying this approximation 
to Px{Sx) in the logarithmic function in the last term in 
equation (113, we obtain the expression of the variational 
free energy as 

X 

Q?}] := E + E 

i^V a—1 

(12) 

P i^V a=l 

where the functionals, Vi[Qi] and Vint[Q“], are defined as 

V^m := - E e^{Si)QAS^) + ^ E 

Si ^ Si 

(13) 

Vi„t[Q“] := E i?mt(5“)Q“(5) + ^ E lnQ“(S“). 

(14) 


ni^vUi=iQ?is^) 


(11) 


In the context of the CVM, the x-replicated free en¬ 
ergy in equation (US is approximated by the minimum 
of the variational free energy in equation (1121) with re¬ 
spect to the marginal distributions {Qi,Q°‘,Qf}, i.e., 
Fx w min{Q._Qc_Q=} Note that, 

at the minimum point, the normalization constraints for 
and the marginal constraints, 

Qfisr)= Y ( 15 ) 


and 


Q?( 5 r)= E ^?“(^“)> ( 16 ) 


should hold. 

In order to minimize the variational free energy with 
respect to {Qi(S'i)}, by using the Lagrange multipliers, 
we perform the variational minimization of 

P{{Q^}] :=E12*[^?*]-E«*(EQ*(^*)-i) 

iev iev Si 
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Hint(S") 




FIG. 3. (a) Factor graph representation of Px{Sx) when n = 2. (b) Cluster decomposition for (a). In the decomposition, the 
three different types of clusters are employed: the clusters consist of S^, of Si, and of 


X 

E QiiSi)-Qns?)) 

i6Va=l Sf Si\{Sf} 

with respect to From the result of this mini¬ 

mization, we obtain 

X 

Q,{S,) oc exp/?(e .(50 + ^ A“( 5 r)). ( 17 ) 

OL — 1 

The Lagrange multipliers, = {A“(5'i) | a = 

1, 2 ,..., a:}, are determined such that they satisfy equa¬ 
tion m- By substituting equation (13 into equation 
(US, while noting the marginal constraints in equation 
(HU, we obtain the partially minimized variational free 
energy, 

_;,RCVM[{Qa^ga}] min ga 

I J 

as 

X 

g?}] = E extr{ ^ ^ Af{Sr)Qf{Sr) 

i^v ' a=l Sf 

- En f dhp^{h)f['^expP{(j)i{S^,h) 

^ a=l Sf 

X ^ X 

+Afisr ))} + ^ Vintm - a E E( 18 ) 

Q—1 i^V a—1 

where the notation “extr” denotes the extremum with 
respect to the assigned parameters. 

B. Bethe Approximation and Replica Symmetric 
Ansatz 

The functional Vint[g“] can be interpreted as the vari¬ 
ational free energy for the interaction term of the original 
system (see equation (HI)- Since this variational free en¬ 
ergy is intractable in general, we approximate it by the 


Bethe approximation. As in equation ([31), we approxi¬ 
mate g“(S“) by 

(n.ev Qfisr)) ( Qusr,s-)) 

Uii,y^EQns?)Qnsf) 

(19) 

where the distributions, Q°^{Sf) and Qfj{S°',S^), are 
the marginal distributions of Q°‘{S°‘), so that the 
marginal constraints, 

Y,QUsr,sf) = Qf{sf) 

Sf 

and 

j2QUsr,sf) = Qnsr), 

Sf 

are satisfied. By applying equation (fT 3 ]) to g“( 5 '“) in the 
logarithmic function in the last term of equation da , we 
obtain the Bethe approximation of Vint[Q'^] as 

^=-E E 

Kilessf.sf 

+ + ^ E (^ 2 [g“,]-Hi[g“]-Hi[g“]). 

^ ieV ^ 

( 20 ) 

By substituting the Bethe approximation in equa¬ 
tion (EOl) into equation (113, we obtain the Bethe ap¬ 
proximation of equation (fT^ : ^Qf}] ~ 

After the Bethe approximation we 
make the replica symmetric (RS) assumption d, [T 4 I in 
■^x^^[{gf I gi^j}]) subsequently, by taking the limit 

as a: —>■ 0 , we finally reach the variational free energy 
expressed as 

= ^extr{ ^ A,(5,)g*(^,) 

iev ' Si 
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-iy Ahpi{h) \ny^ exp 13 {(j)i{Si,h) + A,(5*,))| 

- E E 

+ 4 E (^2[Q,.,]-Hi[Q.]-Hi[Q,]). (21) 


The detailed derivation of this variational free energy is 
shown in appendix We expect that the minimum of 
this variational free energy is the approximation of the 
quenched average of the Bethe free energy in equation 
( 0 ): [Fhetheih, I3)]h « where 


pLBP(RS) , 
min 


min 




( 22 ) 


At the minimum of the variational free energy, 
jrLBP(RS) normalization constraints, 


= E = 1 , (23) 

Si Si,Sj 

and the marginal constraints, 

Y,Q^ASi,S,)=Q,iS,) (24) 

Si 


and 

Sj) = Qi{Si), (25) 

Si 


should hold. 


C. Message-passing Equation 


the two-vertex marginal distributions {Qij{Si, Sj)} are 
obtained as 


Qij{S„Sj) oc Qi{S,)Qj{Sj) exp (l3tfj,j{S„ Sj)) 

X (29) 

The detailed derivation of equations (P5)) "- (l^^ is shown 
in appendix 11 The order of the computational cost of 
the proposed message-passing equation is equivalent to 
that of standard LBP in section IIIBI 

When the distributions of the random fields are Dirac 
delta functions, Pi{h) = 5{h — hi), namely, the fields are 
not the quenched parameters but the fixed parameters, 
the method presented in equations (l^ - (l^ is reduced 
to the standard LBP in equations (IS])-©. In this case, 
in equation (l22l) is equivalent to the Bethe free 
energy Fbethe{h,A- 

After numerically solving the simultaneous equations 
in equations (IMJ-dMI), by substituting the solutions, 
{Qij{Si, Sj)}, and {Ai(S'i)}, into equation 
m, we obtain the minimum values of the variational 
free energy, and regard it as the approxima¬ 

tion of [Pbethe(^,/3)]h. in equation ([5]). 

For the moment, we suppose that the function 
4>i{Si,hi) can be divided as 4>i{Si,hi) = (j)f^\Si,hi) + 
(j)[^\Si). The variations in the quenched average of the 
Bethe free energy in equation ([5]) with respect to 4>[^\Si) 
and 'ipij{Si, Sj) are 


^[Abethejh 


-[HSi)]h 


(30) 


and 


^[Abethe]h. 

S'tpij{Si,Sj) 


[b^As^,SJ)]h, 


(31) 


In this section, we show the message-passing equation 
for minimizing the variational free energy in equation 
m obtained in the previous section. 

The message-passing equation for our method is ob¬ 
tained as 


respectively, which are the quenched average of the beliefs 
obtained from LBP. On the other hand, the variations in 
with respect to and Sj) are 

obtained as 


Si 

(26) 

<3.(S.) = f 

J + A{Si)) 

(27) 

l3Ai{Si) = ln/ife^.(5.). (28) 


kGdi 


The quantity p,i^j{Sj) is the message from vertex i to 
vertex j, and the Lagrange multipliers {Ai(S'i)} in equa¬ 
tion (I28p satisfy the extremal conditions in the first term 
in equation (EH). By using the messages and {(5i(<S'i)}, 


cpLBP(RS) 




-QiiSi) 


(32) 


and 

rpLBP(RS) 

= (33) 

respectively. By comparing equations (1501) and (lOTl) with 
equations El and (l33l) . it can be expected that, if 
is a good approximation of [Fbethe(/i,/3)]^, the 
marginal distributions, Qi(Si) and Qij{Si, Sj), are also 
good approximations of the quenched averages of the be¬ 
liefs, [ 6 i(S'i)]^ and [bij{Si, Sj)]h, respectively. 
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D. Numerical Experiment 

In this section, we describe the evaluation of the va¬ 
lidity of our method by using numerical experiments. In 
the experiments, we used the model expressed as 

P{S\h) = ^e^p[Y.^^^^+ E 

i&V {ijyeE 

(34) 

which is defined on a certain graph, where Si takes q real 
values in the interval [—1,1] as Si € {2S/{q — 1) — 1 | 
S = 0,1,2,..., q - 1}, i.e.. Si S {-1-1, -1} when q = 2, 
Si e (-1-1,0, —1} when q = 3, and so on. The fields h are 
i.i.d. random fields drawn from the Gaussian distribution 
with a mean of zero and a variance of a^, Nihi \ 0, a^). 

We compared the free energy per variable obtained by 
our method, ^ jn, with the quenched average of 
the Bethe free energy (per variable) shown in equation 
®, which was obtained by numerically averaging the 
Bethe free energy in equation (j?]) over the random fields, 
and we compared the behaviors of the quenched average 
of the magnetizations obtained by the two different meth¬ 
ods. The magnetizations obtained from LBP and our 
method are given by Mlbp := n~^ J2ievJ2si >S'i[^i(<S'i)]^ 
and Mrlbp := n~^ J2i(^vJ2si ^iQiiSi), respectively. 

1. Square Lattice 

We show the results when the model in equation (IMl) 
is defined on a graph of an 8 x 8 square lattice with 
the free boundary condition and when all of the interac¬ 
tions are unique, Jij = J. Figures HHS] show the plots 
for g = 2, q = 3, and q = 4. “LBP” represents the 
results obtained by the numerically averaged Bethe free 
energy, and “RLBP” represents the results obtained by 
our method. Each plot of LBP is numerically averaged 
over 10000 realizations of the random fields. In al¬ 
most all cases, the results of our method are consistent 
with the numerically averaged Bethe free energies, as ex¬ 
pected. However, in the cases of large J, mismatches 
between the two methods are observed. 

Figure [7] shows the plot of the quenched average of 
the magnetizations, Mlbp and Mrlbp, when the model 
shown in equation (1341) is defined on a graph of a 14 x 14 
square lattice with periodic boundary conditions and 
when q = 2. We observe that the two methods show 
the different nature of the magnetizations. The magneti¬ 
zation obtained by standard LBP continuously increases 
with the increase in J, whereas that obtained by the pro¬ 
posed method drastically increases, like a first-order tran¬ 
sition, around J « 0.88. This different physical picture 
probably causes the mismatches between the two meth¬ 
ods in the cases of large J in figures HHS] 

Our formulation allows for the approximate evaluation 
of the quenched average of Bethe free energy over the 
random fields with disordered interactions. We show the 


results when the model in equation (1341) is defined on a 
graph of a 14 X14 square lattice with free boundary condi¬ 
tions and when the interactions {Jij} are independently 
drawn from A/"(J^- | 0,(5^). Figure |S] shows the results of 
the free energies versus ct for q = 5. Each plot obtained 
by LBP is numerically averaged over 100 realizations of 
the random fields and over 200 realizations of the interac¬ 
tions, and that obtained by our method is averaged over 
200 realizations of the interactions. Since the error bars 
of our method are quite small compared to LBP, we omit 
them in the figure. The results obtained by our method 
are consistent with the numerically averaged Bethe free 
energies. 

To see the effect of the disorder in the interactions 
on the behavior of the magnetization, we show the plot 
of the quenched average of the magnetizations when the 
model shown in equation (1341) is defined on a graph of 
a 14 X 14 square lattice with periodic boundary condi¬ 
tions and when q = 2 and the interactions {Jij} are in¬ 
dependently drawn from N{Jij \ c, d^) in figure O We 
observe that the magnetizations obtained by our method 
show the first-order transition, as in figure [71 However, 
the values of the magnetizations after the transition are 
quite small compared to those in figure [T] 


2. Random Regular Graph 

A random regular graph (RRG) is a random graph in 
which the degrees of all vertices are fixed by the con¬ 
stant d. Figure [TUI shows the results when the model in 
equation (|M)) is defined on an RRG with 200 vertices 
and d = 4 and when Jij = J and q = 2. Each plot 
obtained by LBP is numerically averaged over 100 real¬ 
izations of the random fields and over 100 realizations of 
the structure of graph, and that obtained by our method 
is averaged over 100 realizations of the structure of graph. 
Since the error bars of our method are quite small com¬ 
pared to LBP, we omit them in the figure, as in figure 
[51 The behaviors of the quenched magnetizations, Mlbp 
and Mrlbp, obtained by the two methods in this case 
are shown in figure 1111 The behaviors of the quenched 
magnetizations in this figure are similar to those shown 
in figure [T] 

LBP is asymptotically justified on an RRG Q, be¬ 
cause an RRG is quite sparse. Therefore, we can expect 
that the results obtained by LBP are close to the exact 
solutions. Except for the RS assumption, our method 
consists of two approximations: the approximation in 
equation (HU and the approximation in equation (HU. 
Since the latter approximation is the Bethe approxima¬ 
tion, it can be justified on a sparse graph such as an RRG. 
This suggests that the mismatch between the two meth¬ 
ods in the right panel in figure |T0] is mainly caused by 
the first approximation, and that the first approximation 
produces the metastable state that causes the first-order 
transition in figure ITT] 

As in section IHID H we again see the case with the 






FIG. 4. Quenched Bethe free energies per variable for q = 2. The left panel shows the free energies versus cr with J — 0.2, and 
the right panel shows the free energies versus J with a = 1. The error bars are the standard deviation. 




FIG. 5. Quenched Bethe free energies per variable for q = 3. The left panel shows the free energies versus cr with J — 0.2, and 
the right panel shows the free energies versus J with cr = 1. The error bars are the standard deviation. 


disordered interactions. Figured shows the plots of the 
quenched Bethe free energies versus cr for g = 5 when the 
model in equation (l34l) is defined on an RRG with 200 
vertices and li = 4 and when the interactions { Jij} are in¬ 
dependently drawn from \ 0, 5^). Each plot in the 

figure is obtained in the same manner as the case in fig- 
ure[8] As is the case in figurelH the results of our method 
are consistent with the numerically averaged Bethe free 
energies. Figure IT^ shows the plot of the quenched av¬ 
erage of the magnetizations when the model in equation 
(131 is defined on an RRG with 200 vertices and d = A 
and when q = 2 and {Jy } are independently drawn from 
Af{Jij I c, (5^). It can be observed that the transition 
of the magnetization obtained by our method is nearly 
continuous with the increase in the magnitude of the dis¬ 
order. This suggests that the disorder in the interactions 
violates the metastable state of the quenched Bethe free 
energy obtained by our method in the case of an RRG. 


E. Exactly Solvable Case — Ferromagnetic 
Mean-field Model in Random Fields 

In this section, we consider the ferrom^netic mean- 
field model in random fields expressed as |6| 

P{S I h.,/3) (xexp|3(^'^(j){S^,h^) + ^Y^g{S^)g{Sj)^, 

iGV i<j 

(35) 

where the second summation represents the summation 
over all of the distinct pairs of vertices, and {hi} repre¬ 
sents the i.i.d. random fields drawn from the distribu¬ 
tion Phihi). By using the Hubbard-Stratonovich trans¬ 
formation, the free energy (per variable) of this model, 
f ■= F{h,f3)/n, can be expressed as 



+ h,) + mg{S) - } 

” i&V S ” 

2n/3 2 tt 

For a sufficiently large n, we can compute the integration 
by using the saddle point method and the summation 
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FIG. 6. Quenched Bethe free energies per variable for g = 4. The left panel shows the free energies versus cr with J — 0.2, and 
the right panel shows the free energies versus J with a = 1. The error bars are the standard deviation. 



FIG. 7. Quenched magnetizations versus J, where q = 2 and 
(7=1. 


over i £ V in the exponent by using the law of large 
numbers, so that the free energy can be expressed as 


where cq and Ci are constants unrelated to Si. From 
this equation, we ensure that all of the messages are con¬ 
stants unrelated to S for a sufficiently large n, because 
^j^i(Si) = exp(ci-l-0(n“^)) ~ exp(ci). Therefore, from 
equations (051) and (1551) . we obtain 

Ai(50 = ^ '^g{Sj)Qj{Sj) + C 2 = mg{Si) + C 2 

j^di Sj 

(39) 

for a sufficiently large n, where C 2 is a constant unrelated 
to S'i, and we redefine m := n~^ Siev Ssi 9{Si)Qj{Si). 
By substituting equation (l3^ into equation (l27l) , we ob¬ 
tain the same expression for m as equation dSZl). Since 
all of the messages are constants, from equation (051) . we 
have 


f = j dhph{h) In'^exp ^{(j){S, h) + mg (S)) 

(36) 


in the thermodynamic limit, where m is the solution to 
the saddle point equation; 


m = 


,, f, s J2s9{S) exp/3(^(5', h) + mg{S)) 
Es exp/3(0(S',/i)-b 7713(5)) 


(37) 


Since the free energy in equation (1551) does not depend 
on {hi}, the quenched average of the free energy over the 
random helds, [F]h/n, coincides with equation (1561) . 

Since 'ipij{Si, Sj) = g{Si)g{Sj)/n, the message-passing 
equation in equation (051) can be expanded as 


\nfj.j^i{Si) 

= ln^Qj(5j)exp (l3'ipij{Si, Sj))fii^j{Sj)~^ + cq 
Si 

P9{S^)Esi9iSi)QiiS,)^^,^,iS,)-^ 

n +ci+0[n ), 

(38) 


By substituting this equation and equation dM]) into 
equation (OT1) . we find that the equaiity / = 
holds in the thermodynamic limit. From this result, we 
can conclude that our message-passing method can ex¬ 
actly compute the quenched average of the free energy of 
the model in equation (I35|) in the thermodynamic limit. 


IV. APPLICATION TO BAYESIAN IMAGE 
RESTORATION 

In this section, we describe the estimation of the sta¬ 
tistical performance of the Bayesian image restoration 
system using LBP by using the proposed method. In 
images, each pixel is allocated in a two-dimensional grid 
and has an intensity value corresponding to the color at 
its position. For an image J G {/i | i G P}, the entry 
h G {0,1 ,..., <7 — 1} represents the intensity of the ith 
pixel. 

Suppose that the original image I is degraded to h 
through a specific noise process. In the Bayesian image 
restoration system, we assume that the original image is 
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FIG. 8. Quenched Bethe free energies per variable versus a for q = 5. The left panel shows the free energies when 5 = 0.2, and 
the right panel shows the free energies when <5 = 0.4. The error bars are the standard deviation. 



FIG. 9. Quenched magnetizations versus c for 5 — Q.2 and 
5 = 0.5 when g = 2 and (7 = 1. 

the sample drawn from the specific distribution Ppri(S) 
known as the prior distribution. We observe only the de¬ 
graded image as the input, and we want to estimate the 
original image from the input degraded image. To do the 
image restoration, we compute the posterior distribution 
of the original image, Ppost(S^ | h) oc Piike(^ | S)Ppri{S), 
where Pii\ie{h \ S) is the noise process referred to as 
the likelihood. We use the posterior distribution to pro¬ 
duce the restored image as the output. The scheme of 
Bayesian image restoration is shown in figure 1141 

A. Bayesian Image Restoration nsing LBP 

For the original image /, suppose that the degraded 
image h is generated by adding additive white Gaussian 
noise, i.e., hi = li + rji, where pi is the random noise 
drawn from the Gaussian distribution \ 0, tr^). 

We define the prior distribution of the images S as 

Ppri(5 I a) oc exp ^ ^{Si,Sj)^ , (40) 

{hlleB 

which is defined on a square lattice corresponding to the 
configuration of pixels. The energy function ^(Si, Sj) de¬ 


fines the relationship among neighboring pixels, and it of¬ 
ten takes a form that emphasizes the smoothness among 
neighboring pixels, e.g., ^{Si,Sj) = — |5'i — 5'i|. The posi¬ 
tive parameter a controls the strength of the smoothness. 
For the original image, the distribution of the degraded 
image, i.e., the likelihood, is expressed as 

Piike(/i \S = I,a‘^) = l[M{h, I h,a^). (41) 

iGV 

From equations (1401) and dm, given a specific degraded 
image h, the posterior distribution of the original image 
is obtained by 

Tpost {S I h,a,a'^) (xPiike{h \ S,a'^)Ppri{S \ a) 

ocexp(-^ -ba ^ S,{Si,Sj)\. (42) 

iey ^ 

The posterior distribution is the special case of equation 
(ED- The degraded image is regarded as the random fields 
in the posterior distribution. 

In maximum posterior marginal (MPM) estimation, 
the zth pixel of the restored image is obtained by 

5'jMPM := arginaxPpost(>S'i \h,a,a‘^), (43) 

where Ppost(<5'i | h,a,a‘^) is the marginal distribution of 
the posterior distribution in equation dni) El. In prac¬ 
tice, we approximate the true marginal distributions by 
the beliefs obtained by LBP for the posterior distribu¬ 
tion, 

^MPM ^ ^MPM(LBP) 

where the belief hi{Si) is obtained by LBP described 
in section IIIBI with (j)i{Si,hi) = —{Si — hi)‘^f2a‘^, 
ipij{Si, Sj) = a^{Si,Sj), and /3 = 1. The performance 
of the restoration is often measured by the mean square 
error (MSE) between the original image and the restored 
image, 

D{I,h,a,a^) := (45) 

^ iev 
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FIG. 10. Quenched Bethe free energies per variable for q = 2 on the RRG. The left panel shows the free energies versus a with 
J = 0.2, and the right panel shows the free energies versus J with ct = 1. The error bars are the standard deviation. 



FIG. 11. Quenched magnetizations versus J on the RRG, 
where q — 2 and cr = 1. 


This MSE is for the specific input degraded image that is 
drawn from the likelihood in equation (14111 . and it should 
take different values for different degraded images in gen¬ 
eral. 

For the original image /, we attempt to estimate the 
average value of the MSE in equation (1^ over all pos¬ 
sible degraded images, 

Zlav(i', a, o-^) := [D{I, h, a, cr^)] (46) 

where Pi{h) = Af{h \ Ii,a^), because the average value 
corresponds to the statistical performance of the pre¬ 
sented Bayesian restoration system for the original image. 
By considering the input degraded image as the quenched 
parameter, the right-hand side of equation (H51) can be re¬ 
garded as the quenched average of the MSE obtained by 
LBP. Thus, we approximate it using the message-passing 
method presented in section UlI Cl 

Using the proposed message-passing method, we ap¬ 
proximate equation (HSI) as 

L>av(-f,a,cr^) ~ “ X! I /i,cr^)(/i - 

” i&V •' 

( 47 ) 


where 

ri(h) = a.rgmax + Ai{S)). 

The detailed derivation of this approximation is shown 
in appendix [Cl 


B. Numerical Experiment 

In this section, we describe the estimation of the per¬ 
formance of Bayesian image restoration for the 64 x 64 
original colored images shown in figure 1151 Colored im¬ 
ages consist of three different channels: red, green, and 
blue (RGB) channels, I = {Jred, ^green, ^biue}- The pix¬ 
els in each channel in the original images have eight 
intensities, i.e., q = 8. In the following experiments, 
we assume that the three different channels are inde¬ 
pendently degraded by the same noise process in equa¬ 
tion dSD, and we restore the generated degraded im¬ 
ages, h = {Ked,hgreen,huue}, by Separately applying 
the Bayesian image restoration based on the posterior 
distribution in equation to the RGB channels. The 
value of cr^ used in the restoration is the same as that 
used in the noise process, and we use the same values 
of the parameters, a and cr^, in the restorations for the 
three different channels. The total MSE of the restora¬ 
tion is obtained by the average of the MSEs of the RGB 
channels, i.e., D{I,h,a,a‘^) = {i4(Jred,/^red, a, cr^) -|- 

d(Egi.een5 ^greent (Xj fJ ) -(- ZI(Edluei ^bluei ^ )}/^' 

For the original colored images in figure 1151 we eval¬ 
uated the average of the MSE, Uav(J, a, cr^), by using 
two types of methods: LBP and the proposed analyti¬ 
cal method in equation (1471) . In LBP, we approximated 
Uav(J, a, cr^) by the sample average of D{I,h,a,a^) 
over 10000 different degraded images, which are gener¬ 
ated from the original image I through the noise process 
shown in equation m- 

At first, we show the results obtained by setting the 
function as ^{Si,Sj) = —{Si — Sj)‘^/2 in the 

prior distribution in equation (HiH) . In figures (TH] and [T71 
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FIG. 12. Quenched Bethe free energies per variable versus cr for g = 5 on the RRG. The left panel shows the free energies when 
S — 0.2, and the right panel shows the free energies when S — 0.4. The error bars are the standard deviation. 




FIG. 15. Original 3-bit colored images I: (a) lenna, (b) par¬ 
rots, and (c) sailboat. 


FIG. 13. Quenched magnetizations versus c on the RRG for 
S = 0.2 and 5 = 0.5 when g = 2 and a = 1. 


PpAS) 


Plike{h I S) 

additive noise 



estimate 


post 


is\h) 


posterior 


FIG. 14. Scheme of Bayesian image restoration. 


we show the plots of Dav{I,ci,o''^) versus a and a, re¬ 
spectively, obtained by the two different methods: the 
sample average of the LBP restorations, “LBP,” and the 
proposed analytical method, “RLBP.” Figure [16] shows 
the plots of cr, cr^) versus a when tr = 0.5, and 

figure [Tfl shows the plots of a, cr^) versus a when 

a = 0.4. Each plot obtained by LBP restoration is aver¬ 
aged over 10000 realizations of the stochastically gener¬ 
ated degraded image. We can observe that the results 
obtained by our method in equation (H71) are in good 
agreement with those obtained by LBP restoration. 

Next, we show the results obtained by setting the func¬ 
tion ^j) as ^{Si,Sj) = —15—Figure [IHI shows 


the plots of Dav{I, a, cr^) versus a when a = 0.5, and fig¬ 
ure m shows the plots of Dav{I,C(,a^) versus a when 
a = 0.4. The results obtained by our method are in 
good agreement with those obtained by LBP restoration 
in figures Uni and |T71 


V. CONCLUSION AND REMARKS 

In this paper, we proposed an analytical method based 
on the idea of the RCVM to approximately evaluate the 
quenched average of the Bethe free energy, obtained by 
LBP, of the pair-wise MRF in equation m in random 
fields h. Since our message-passing-type formulation al¬ 
lows any form of the functions 4>i{Si, hi) and ipij{Si, Sj) 
in the Hamiltonian and allows any form of the distri¬ 
butions of the random fields {piihi)}, except for the 
cases where correlations among fields exist, the proposed 
method is applicable to a wide range of applications in 
physics and in computer science. In the argument in sec¬ 
tion IIIIEl we found that this approximation is justified 
in the ferromagnetic mean-field model in random fields. 

Although the results obtained by our analytical 
method in almost all cases were consistent with those 
obtained by the numerical method, as seen in the artifi¬ 
cial model presented in section llll PI and in the Bayesian 
image restoration presented in section lIVBi some mis¬ 
matches, especially in the behaviors of the phase transi- 
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FIG. 16. Plots of Da,v{I,a,a^) versus a when a = 0.5 and ^{Si,Sj) = —{Si — Sj)^/2 for the original images shown in figure 
[13 the left, middle, and right plots show the results for figures [I3;a),[I3b), and (He), respectively. The error bars are the 
standard deviation. 




FIG. 17. Plots of Dav(/, a, (T^) versus cr when a = 0.4 and ^{Si,Sj) = —{Si — Sj)^/2 for the original images shown in figure 
[13 the left, middle and right plots show the results for figures [I3;a), [I3b), and [He), respectively. The error bars are the 
standard deviation. 


tions of the magnetizations in section IIII Di between the 
results obtained by our method and by the numerical 
method were observed. As mentioned in section [III D 21 
the mismatches are considered to be mainly caused by 
the approximation based on the CVM in equation (EH- 
However, a detailed understanding of its mathematical 
meaning is still unclear, even though it is one of the most 
important components of our method, because the ap¬ 
proximation is for the system replicated by the replica 
method; thus, the development of a mathematical un¬ 
derstanding of how it affects the present method is not 
straightforward. It should be considered in future stud¬ 
ies. 


We employed the Bethe approximation in equation 
m for the purpose of evaluating the quenched average 
of the Bethe free energy. Our approximate framework, 
however, allows us to employ other mean-field approx¬ 
imations instead of the Bethe approximation, and the 
resulting form can be expected to be an approximation 
of an employed approximation method. Thus, in accor¬ 
dance with the presented framework, it is expected that 
we can evaluate the quenched averages of the employed 
mean-field methods for random fields. 


Appendix A: Derivation of Equation (1211) 

By substituting equation ([20l) into equation (HU) , we 
obtain the Bethe approximation of Q“}] as 

X 

Qtj}] ■■= E ^ ^ ^ns-)Qns-) 
--In / +AnsT))} 

+ E Qi,}] - 4 E E (Ai) 

a—1 ze V a—1 

At the minimum point of this variational free en¬ 
ergy, we make the RS assumption; that is, the rela¬ 
tions Qf{Si) = Qi{Si) and Qfj{Si,Sj) = Q,j{Si,Sj) 
hold for any a. Under this assumption, the minimum of 
(5“j}] is equivalent to the minimum of the RS 
variational free energy expressed as 

X 

^extr{ E E 

igy * a=l Sf 
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FIG. 18. Plots of Dg.v{I, a, a^) versus a when a = 0.5 and ^{Si, Sj) = —\Si — 5'j| for the original images shown in figure fTSl the 
left, middle, and right plots show the results for figures [TSJ a), [T^b), and ll5f ci. respectively. The error bars are the standard 
deviation. 



FIG. 19. Plots of Db.v{I, ce, cr^) versus a when a = 0.4 and ^{Si, Sj) = —|Si — Sj \ for the original images shown in figure fT^ the 
left, middle, and right plots show the results for figure [T^ a'l. IlSl bi. and IlSf ci. respectively. The error bars are the standard 
deviation. 


-4ln / d/ipi(/i) ^exp/3(0i(5'“,/i) + A“(5'“))| 

^ '' a=l S“ 

+ (A2) 

^ i&V 


From the convexity of the first term of this equation with 
respect to Aj and the extremal conditions of Aj, it can 
be ensured that the relation Af{Si) = Ai{Si) holds for 
any a at the extremal points of Aj. Therefore, equation 
(IA2I) can be reduced to 


and we regard the minimum of this variational free energy 
as the Bethe approximation of the true cc-replicated free 
energy in equation (HQl). 

From equations ( 0 , cni), and (ESI, our desired varia¬ 
tional free energy is obtained by 


P X 


(A4) 


This leads to equation (1^ . 


jrLBP(RS)j|Q.^g_|j = ^extrjx^ Ai(S'0Qi(S'i) 

iev ‘ Si 

-iln j dhpi{h) exp (^xln'^exp P{(j)i{Si,h) + Ai{Si))^'^ 

P J Si 

+ - I ^ Hi[Q,], (A3) 

P iev I 


Appendix B: Derivation of the Message-passing 
Equation 

To perform the conditional minimization of the varia¬ 
tional free energy in equation (EH, we use the Lagrange 
multipliers as 


^LBP(RS) _^LBP(RS) ^ a. ( ^ Q^S^) - l) 

iev Si 


Si.Sj 
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- E 

b,j}&E 


{ E (E + E ^^ASj){ E Q^AS^, s,) - Q,(S ,))}, 


where the Lagrange multipliers, {oi, bij} and {MjiSj}}, 
correspond to the normalization constraints in equa¬ 
tion (|23l) and the marginal constraints in equations 
(I24|) and (|25ll . respectively. From the minimum condi¬ 
tions of Qij}] with respect to Qi(Si) and 

Qij{Si, Sj), we obtain 

Q,{S,) oc exp E (a,(, 5,) + ^ Afe,,(5,)) (Bl) 

I kGdi 

and 

Qi,j{Si,Sj) oc expfi{tp^j{Si,Sj) + XijiSj) + Xj^i{S^)), 

(B2) 

respectively, where the notation | • • • | denotes the num¬ 
ber of entries of the assigned set. By introducing the 
messages in the form 

pLj^i{Si) := exp I ^ (^Ai(5'i) - /3“^ ln(5^(S'i) 

+ ^ Afc,,(5,)) -/3A,„(5,)}, 

kGdi 

we obtain 

exp (pXj^iiSi)) = (5i(5'i)exp ( - (3Ai{Si)) 

X Y[ (B3) 

feG9i\{j} 

From equations (ED and (IB3|) . we obtain the relation 
PAi{Si) = ^ hifik^i{S^) + a, (B4) 

kGdi 

where Ci is a constant unrelated to Si. Since the value of 
Ci does not affect our results, without loss of generality, 
we set {ci} to zeros. 

By substituting equations (IBlI) and (IB2I1 into the 
marginal constraints in equations (IMl) and (1^ and using 
equation (IB3I) . we obtain the message-passing equation 
as 

(X '^Qj{Sj)expP{ - Aj{Sj) + i’ij{S^,Sj)) 

Sj 

X l^k^j{Sj) 

k&dj\{i} 

OC E (*^4: )) Mi—) 

Sj 

(B5) 

Here, from the first line to the second line of this equa¬ 
tion, we use the relation in equation (lB4l) . From the 
extremal conditions for {Ai(S'i)} in the first term in equa¬ 
tion (ED, we obtain equation ED- From equations (IB2I) , 
(IB3|) . and (lB4l) . the two-vertex marginal distributions 
can be expressed as equation (1^ . 


Appendix C: Approximation in Equation (14711 

Using the belief obtained by LBP and any real value 
7 , we define the new distribution as 


C^{Si I 7) := 


bliS^) 


Es.bliS.) 

hi) + Y.]^di 

Esi exp7(/3(^i(S'i, hi) -f Yjjddi lnM,-^i(5',)) 


• (Cl) 


Assuming b'l{Si) has a unique maximum with respect to 
Si, we obtain the equality 


(argmax5i(S'i))^j = lim ^ S'f | 7 )](C2) 

Si \h J-ioo ^"• 


for any k. 

As mentioned in section UlICl the distribution Qi(Si) 
in equation ED, 


QiiSi) = / dhpi{h)qi{Si I h), 


where 


QiiSi I h) := 


exp^{(j)i{Si,h) + Ai{Si)) 

expP{(j)i{Si, h) + Ai{Si)) ’ 


is regarded as the approximation of the quenched average 
of the belief, [&(5'i)]h. According to equation (ICll) . we 
define the distribution as 


Pz{Si I h,7) := 


QiiSi I h)'> 


T,s, I h)7 

exp"fl3{4)i{Si,h) + Ai{Si)) 

J2si exp7/3(Si,/!) +Ai{Si)) 


As mentioned in section UlICl when the distributions of 
the random fields are Dirac delta functions, Pi{h) = 5{h— 
hi), the proposed method is reduced to standard LBP. 
Thus, in this case, since qi{Si \ hi) = bi{Si), the equality 
Pi{Si I hi, 7 ) = Q{Si I 7 ) holds. From this relation, it 
is expected that the quenched average of I l) is 

approximated as 


[m{Si I 7 )]^ - J dhpi{h)pi{Si I h, 7 ). 


Using this approximation, we approximate the right- 
hand side of equation (IC2|) by 


hm | 7 )]^ 


7 —^■OO 
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Ri lim '^Si [ dhpi{h)pi{Si \ h,-f) 

'y—^co I 

Si 

= J dhpi{h)ri{h)^, (C3) 


where 


ri{h) := argmax h) + Ai(5')). 

Using equations (Hll) . (IC2I) . and (IC3I) . we obtain the 
approximation as 


j-^^MPMCLBP) 



is approximated by 

~ - X! f dhpi{h)ri{h) 

^ i&v ^ 

- J dhpi{h)ri{hf'j 

= f dhpi{h){li - ri{h)Y. 


[( 5 : 


,MPM(LBP)\ fc 


)\R.y'dhp,(/r)^(h)^ (C4) 
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